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We investigate how forces spread through frictionless granular packs at the jamming transition. 
Previous work has indicated that such packs are isostatic, and thus obey a null stress law which, 
independent of the packing history, causes rays of stress to propagate away from a point force at 
oblique angles. Prior verifications of the null stress law have used a sequential packing method 
which yields packs with anisotropic packing histories. We create packs without this anisotropy, 
and then later break the symmetry by adding a boundary. Our isotropic packs are very sensitive, 
and their responses to point forces diverge wildly, indicating that they cannot be described by any 
continuum stress model. We stabilize the packs by supplying an additional boundary, which makes 
the response much more regular. The response of the stabilized packs resembles what one would 
expect in a hyperstatic pack, despite the isostatic bulk. The expected stress rays characteristic of null 
stress behavior are not present. This suggests that isostatic packs do not need to obey a null stress 
condition. We argue that the rays may arise instead from more simple geometric considerations, 
such as preferred contact angles between beads. 

PACS numbers: 45.70.-n, 62.20.D-, 05.10.-a 



I. INTRODUCTION 

The study of forces in static granular bead packs has 
a long tradition, owing to the rich panoply of observed 
phenomena. Extensive work has been done to study the 
distribution of of the inter-bead contact force strengths 
PHI] and angles [5] , the presence of load bearing arches 
[5J [7] , and the scaling of various properties such as pres- 
sure [51 [S], effective shear modulus [HI US], and length 
scales [TT1 [T^] . The large assortment of interesting be- 
haviors is due partly because the force response depends 
on the particular microscopic arrangement of the beads. 
Altering the geometric structure of the pack may have 
almost no effect at all, or it may cause huge changes to 
the bulk properties of the material. 

In this paper, we do not attempt to discuss how 
the forces will spread in an individual disordered pack. 
Rather, we attempt to describe the force propagation 
using a continuous, coarse-grained stress tensor by av- 
eraging over many different packs and finding the mean 
response. Such an approach is possible for hyperstatic 
granular packs [El [14]. However, it is not guaranteed to 
work for the isostatic packs we consider here. 

Broadly speaking, the set of possible static granular 
packs can be divided into three types based on their sta- 
bility. Hyperstatic systems are the most common. The 
number of contacts between beads is large, and if an ex- 
ternal force is applied to the pack, forces can be generated 
across these contacts to stabilize the load. Hypostatic 
packs, on the other hand, do not have many contacts be- 
tween beads. They tend to be less dense, and are unable 
to support some external loads. Such packs are charac- 
terized by the existence of so-called floppy modes, which 
are internal rearrangements of the bead positions which 
do not require any energy. Isostatic packs exist at the 
transition between these two types. They are stable in 
the sense that they can support any infinitesimal load 



provided. However, if even one contact is removed, they 
become hypostatic: a floppy mode is created, which typ- 
ically causes a rearrangement of beads across the entire 
pack. 

The stability of a pack can be estimated by comparing 
the number of beads in it with the number of contacts. 
In ci-dimensional granular packs consisting of N beads, 
static packs must satisfy dN force balance equations. 
Without friction, the forces on any bead are either due 
to contacts with its neighbors or are applied externally. 
If the beads have, on average, Z contacts each, then the 
total number of contacts in the pack is n = NZ/2. If this 
is less than dN, then the system is undcrdctcrmincd. As 
discussed above, this happens in hypostatic packs. If 
n > dN, then there are more contacts than necessary. 
This is typically the case in hyperstatic packs. On the 
other hand, if n = dN, then often the pack is isostatic. 
In two dimensions, this means that the average contact 
number is Z — 4. Because an isostatic pack has the same 
number of contacts as force balance equations, one can 
solve uniquely for all of the contact forces once the ex- 
ternal forces acting on the system are known. Note that 
this isostatic contact number is necessary but not suffi- 
cient for the system to actually be isostatic; for a true 
test, one must check for the presence of floppy modes. 

Isostatic packs behave quite differently than their hy- 
perstatic brethren. Our goal is to find a continuous, 
coarse-grained stress tensor a a p at all points in the pack. 
In this continuum limit, the stress tensor should obey the 
continuity equations 

d a <r a fj = Fp (1) 

where F is the external force supplied to the system. 
These equations express static force equilibrium and, in 
two dimensions, produce two equations for the three com- 
ponents of the stress tensor. These two equations are not 
sufficient to determine the continuum stress field; a third 
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piece of information is necessary. In a conventional elas- 
tic medium, this final equation is the constitutive stress- 
strain law, which reduces to the well known Hooke's Law 
for small deformations |T51 Section 4]. This is not sensi- 
ble in an isostatic medium where no contacts need to be 
compressed; as discussed above, on a microscopic level, 
force balance is actually sufficient to determine the re- 
sponse in any particular case. One proposal for a final 
continuum equation in the isostatic case is a linear re- 
lationship between the components of the stress tensor 



For a quick idea of why the stresses should be linearly 
related in isostatic systems, imagine cutting out a small 
section from a large two dimensional isostatic pack. This 
will result in the creation of several floppy modes, due to 
the removal of the contacts on the sides. Now apply a 
stress to two opposite sides of this section. In an elastic 
material, there would be some resistance — energy could 
be stored in the resulting compression, and the sample 
would try to push back. However, the floppy modes allow 
the beads to move without using energy. The positions in 
two directions are restricted by the applied deformation, 
but the system is still free to move out in the other two 
directions. To confine the sample, we must also push 
on the section in these other directions, by an amount 
proportional to the force exerted on the sides. Thus in 
isostatic packs, we should have 
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This argument does not depend on any packing history, 
and should hold for all isostatic packs. 

With the linear stress relationship established, the 
equations for each of the separate components of the 
stress tensor can be decoupled [16]. Taking derivatives 
of the continuity equations ([I]) gives 
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Plugging in the linear relationship from Equation[2]allows 
us to eliminate a yy . This gives 
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Solving for a xx and a yy in a similar manner gives the 
same hyperbolic differential equation. Note that if /i = 0, 
then this becomes 



(a) 




FIG. 1. Cartoons of dyy as a function of x, a distance y below 
the forcing point. In (a), we show the setup: a point force is 
supplied to the pack, and we look at a yy in a region a distance 
y below it. In (b), we see the expected stress in this region 
for both elastic and isostatic media. The elastic solid has a 
large central maximum below the applied force, whereas the 
isostatic pack has two peaks, corresponding to the positions 
of the light cones. 



This is the wave equation, with the y direction serving 
the role of time. Just as a flash of light will produce 
light cones extending out into spacetime along null lines, 
a point source of force will create lines of stress extend- 
ing through the pack at angles determined by the "speed 
of light," c. For this reason, we refer to the linear re- 
lationship as a "null stress" condition. In keeping with 
this analogy, we will refer to the lines of stress as "light 
cones." In packs with left-right symmetry, the [i = con- 
dition should hold. This is the case in all of the packs we 
examine. Even if /i ^ 0, however, we can still rotate into 
a frame where fj, vanishes and recover the wave equation. 
In this case, the light cones will still exist, but will be 
bisected by some direction other than the vertical. 

This hyperbolic behavior is not seen in elastic media 
or hyperstatic granular packs. In those systems, there 
is a single maximum below the location of the forcing 
point, whereas in a hyperbolic system there is a minimum 
with two peaks to either side, as depicted in Figure [T] 
There are some models where similar rays of force are 
expected in hyperstatic packs OET]. In these models, 
there is a crossover region where at a certain depth the 
positional disorder in the system causes the light cones 
to merge into a central maximum following elliptic rather 
than hyperbolic behavior. This is different than what we 
are describing here, and we do not expect such a central 
maximum 22 . In Section . 



IV we will briefly discuss how 



(d 2 y - c 2 d 2 x )a aP = 0. 



(8) 



the noise caused by disorder obeys different statistics in 
isostatic packs. 

To test this null stress theory, we can supply a small 
force to a bead and see if the stress is, on average, concen- 
trated into light cones. Unfortunately, experimental tests 
have so far proven inconclusive |23|I24] , due to various dif- 
ficulties such as creating packs which satisfy the isostatic 
conditions, measuring the forces in the bulk, eliminating 
friction, and getting enough data to provide a meaning- 
ful average. In the absence of real-world verification, we 
turn to simulation. According to Huygens's principle, 
in n + 1 dimensions, the Green's function for the wave 
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equation has a long tail if n is even. This makes signals 
harder to distinguish than they would be if n were odd 
[251 Chapter 6]. Therefore we work in two dimensions, 
where the light cones are easier to see. The beads that we 
use are frictionless, and interact only with their neighbors 
via radial harmonic contact forces. We will allow these 
contact forces to be tensile. In the next section, we will 
argue that the particular choice of force law should not 
have a large effect on the measurements we are concerned 
with. 

We begin by providing an overview of the systems in 
question in Section [TTJ In Section |TO} we describe previ- 
ous simulations which appear to verify the null stress law, 
and highlight a broken symmetry inherent in the con- 
struction of these systems. We attempt to make a more 
general system in Section |IV[ by creating isotropic iso- 
static packs, and note disagreements with the null stress 
law: not only do we see the response fail to produce light 
cones, but we find that there is no continuum stress field 
at all. In Section|Vj we show that it is possible to restore 
continuum behavior, but produce elastic-like responses 
rather than light cones. In Section |VI| we argue that 
the light cones found in previous simulations may have 
been caused by the contact angle distribution of the packs 
they used, rather than a null stress effect. We then offer 
concluding remarks in Section |VII| 

II. BACKGROUND 

We describe here the methods we use to analyze the 
packs after they have been created. We will work in 
the regime of linear response, where the magnitudes of 
the forces we supply are infinitesimal. Since the system 
is jammed, the positions of the beads do not change in 
response to these forces, and they can be balanced via 
contact forces with other beads. For the system as a 
whole, we can write this as 

Mf = F (9) 

where f is a vector containing the values of the contact 
forces for each of the n pairs of touching beads, F is a 
vector of length 2N containing the x and y components 
of the external force on each bead, and M is a 2N x n 
contact matrix. Typically some sort of hard floor will 
support the pack. This floor is in contact with beads, 
but we do not require the net force on it to be zero. The 
floor can thus balance any net force applied to the bulk 
of the pack. 

Since the system is isostatic, M is square, and for any 
given set of body forces one can uniquely solve for the 
network of inter-bead contact forces which keep the sys- 
tem stable. Furthermore, M depends only on geometry; 
once the contact angles are known, nothing else is re- 
quired to solve the problem. Each row of M represents 
either the x or y components of the contact forces exerted 
on one bead. The row has one non-zero entry for each 
contact that bead has. The value of this entry is given 



by the unit vector pointing from the row's bead to the 
other bead it is in contact with, projected onto whatever 
direction the row corresponds to. This gives each column 
of M a simple structure. Say that column a is due to the 
contact between beads i and j, and is the unit vector 
pointing from i to j. The column will have exactly four 
nonzero entries: in the rows corresponding to bead i's 
x and y components, it will have (riij) x and (nij) v , and 
for the rows corresponding to bead j, it will have {nji) x 
and {jiji) y . Furthermore, = — n^, and |riij| = 1, so 
each column depends on only one unique piece of infor- 
mation. Note that if the contact is between a bead and 
the hard floor rather than between two beads, then the 
column will have only two nonzero entries. In two dimen- 
sional isostatic systems, the average number of contacts 
per bead is four, so as the system size grows, only a van- 
ishingly small fraction of the entries are nonzero. M is 
evidently quite sparse. 

If, for some F, there is no f which satisfies Equation |9j 
then the system is not jammed: no set of contact forces 
can prevent at least one bead from experiencing a net 
force and shifting in response. However, as long as our 
system is isostatic, this will not happen. The positions 
of the beads remain fixed, and so the contact angles can 
never change; M is therefore constant. 

Given an F, some of the elements of f may be negative, 
indicating tensile forces. Therefore, this cannot strictly 
represent a hard sphere system where all contacts are 
purely repulsive. However, Tkachenko and Witten [26] 
have developed an "adaptive network" technique which 
forces all of the contacts in an isostatic system to be 
compressive. They do this by making slight changes to 
the packing geometry so that under a compressive load, 
there are no tensile forces. They discovered that the re- 
sponse in their system followed the null stress prediction 
both before and after this procedure, indicating that the 
existence of light cones does not depend on any restric- 
tions, or lack of restrictions, on the sign of the contact 
forces. Thus for simplicity, we simply allow tensile forces 
in our simulations. There are many examples of physi- 
cal systems where this is a realistic model. For example, 
if neighboring beads are slightly wet, they can form liq- 
uid capillary bridges which hold them together [27], and 
colloids often have attractive interactions due to, for ex- 
ample, depletion [2S]. Often such systems have a yield 
point, where a large enough tensile force will split the 
beads apart. However, in the linear response regime ap- 
propriate for the small forces in question, the contact 
force magnitude should be small enough that this is not 
a concern. 

Instead of using the contact matrix M described above, 
the linear response of the system can also be measured 
using a dynamical matrix T> [5J . While the contact matrix 
depends only on the angles between beads, the dynamical 
matrix assumes a spring constant k between each of the 
connected pairs. With this spring constant, the energy 
of the system is quadratic in the small displacements u 
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each bead has from its equilibrium position: 

£ = k u T Vu 

2 

This T> matrix must be symmetric, and it follows that 
the force the springs exert on each of the beads is —kT>u. 
Thus in equilibrium we can relate the body force on each 
bead to its displacement via F = kVu. This is an equiv- 
alent way of finding the response. To translate between 
the two, note that the spring energy can also be written 
in terms of the contact forces. Then 



£= 2l fTf= 2l( FT(MTrl )( M " F ) 
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Since T> is symmetric and this must be true for all u, 
comparison with Equation [lO] shows that the dynamical 
matrix can be written down solely in terms of the con- 
tact angles: T> = MM T . However, whereas the contact 
matrix only depended on these angles and not the func- 
tional form of the interaction between beads, by using 
the dynamical matrix we are implicitly assuming a har- 
monic interaction. This choice of a linear spring force 
between contacts is not realistic for some physical sys- 
tems. The contact between smooth, solid bodies is more 
accurately described by Hertz's Law in which the the 
force varies as the 3/2 power of the compression rather 
than linearly. However, the properties we study should 
depend only on the isostatic nature of the contact net- 
work, and not on the force law of the contacts. Previous 
simulations have verified that the density of states and 
other features at the isostatic point are unchanged if one 
uses Hertzian rather than harmonic contact forces [TT] , 
Indeed, all of our systems have zero pressure, so any par- 
ticular choice of overlap potential has no chance to cause 
an effect. Once the final bead configuration is reached, 
all analysis for the response can be done in terms of the 
contact matrix AI, which depends only on the angles be- 
tween beads in contact, and not any displacements due 
to a particular force l aw. 

Note that Equation 11 required the use of M^ 1 , which 
does not exist if the system is not isostatic. However, 
even if M is not invertible, it remains true that T> — 
MM T . The proof is more tedious in this case, and as 
such has been relegated to Appendix [A"| 

If each of the beads has some mass m, then mil + 
kVu = 0. Thus the spectrum of V gives the vibrational 
modes of the system: the squared normal more frequen- 
cies ivf are given by the eigenvalues of kV/m. We can use 
these frequencies to measure the density of states, D(uj), 
defined as 



D(u>) = lim 



of LJi such that uj < LUi < lu + Aw) 
2iVAw 



If T> is singular, then some u>i will be 0. Their correspond- 
ing eigenvectors are the floppy modes of the system. 
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FIG. 2. The density of states for an isostatic pack has a 
plateau at low frequencies. A hyperstatic pack, on the other 
hand, has a density of states which goes to zero at small 
frequencies. In this it resembles an elastic solid, where D(uj) ~ 
oj d_1 , with d the spatial dimension. The isostatic pack here 
was created using the method described in Section [IV A| and 
the hyperstatic pack was found by treating all of the nearest- 
neighbor beads in this same pack as being in contact. This 
gives it an average of Z = 5.7 contacts per bead. 



In d-dimensional hyperstatic packs, D(u>) 



,d-l 



for 



small uj. However, the density of states of an isostatic 
pack has a large nonzero plateau at low frequencies, as 
shown in Figure[2j If the pack is only slightly hyperstatic, 
in the sense that the contact number is only a little larger 
than the isostatic value of four, then the density of states 
has been shown to grow linearly from zero until some 
uj* , above which it will exhibit the plateau expected in 
an isostatic pack [TTJ . An isostatic pack will behave 
differently from the hyperstatic one below this uj* . It is 
these low frequencies which should be most important for 
static force propagation. 

In a hyperstatic system M may not be square, so it 
cannot be inverted to give f. To find the equilibrium 
forces in this case, we treat each contact like a harmonic 
spring, as discussed above. The contact force vector is 
then the f which satisfies M i = F while minimizing the 
energy ^f T f from Equation [IT] 

After the contact forces are known, the stress can be 
computed using an established local coarse graining pro- 
cedure O [T3] , which sums up the contributions from all 
pairs of contacting beads: 



<y a fi (x) = fij Rij Kj )a(THj)p 



J ds$(x - Xj + s 



(12) 

Here Ry is the vector pointing from bead i to bead j, 
whose unit vector is ny. Xj is the position of bead j, 
and $ is a non-negative coarse-graining function centered 
at the origin whose integral over all space is 1. Often 
$ is chosen to be a Gaussian or a smoothed Heaviside 
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function. For simplicity, we will use 
$(r) : 



1 

TTW 2 

0; 



r < w 
otherwise 



(13) 



where the width w will often be chosen to be on the order 
of a bead diameter. 

We would like to measure the response in a way that 
takes into account the direction and magnitude of the 
source force, F. The vector (F-cr) has to be proportional 
to the strength of the applied force. Using Equation [T2] 
the constant of proportionality is 

F |fy X) = E (/ ds$ ( x - ^ + (14) 



where 



fi 



(15) 



The response function G contains the information which 
gives the relationship between the contact forces and the 
applied external force. It can easily be computed at 
the location of any contact without reference to coarse- 
graining integrals. This response is similar to the one 
used by Head, Tkachenko, and Witten [3D], but it differs 
in that it is more directly related to the stress. This G 
has units of length, given by the factor of Rij which tells 
the distance between the beads in question. In the fol- 
lowing, we will always report G in units of the average 
diameter of the beads in the pack. 

Many of our results will be presented in terms of G. It 
must be averaged over several different systems to find a 
continuum response. Each instance is individually quite 
noisy, but across several realizations the average will be 
smoother, as shown in Figure [3] 

As a concrete example of how this works in practice, 
consider the small example pack depicted in Figure [4] 
We will write M so that the x and y components of force 
balance for bead i correspond to rows 2i — 1 and 2i re- 
spectively. The columns will represent contacts A-F, in 
order. Then we have 
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To supply a unit force down on bead 1, we use F = 
(0, — 1, 0, 0, 0, 0) T . Since the system is isostatic, we can 
simply solve Equation [9] to get the contact force vector, 
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If the system had been hyperstatic, then we would in- 
stead solve the linear least squares problem |32j of mini- 
mizing |f| 2 subject to the constraints Mi — F. Routines 




FIG. 3. The response, computed using the techniques de- 
scribed in Section III] to a point force (a) in a single pack with 
512 beads and (bjaveraged over 2970 different realizations. 
The packs were generated using the sequential method of Sec- 
tion [III] and for each different pack, several realizations were 
generated by applying the point force to different beads near 
the top of the pack. The response for each realization was 
measured along a line five bead-diameters below that point. 
We choose five bead-diameters here, and elsewhere in this pa- 
per, but the light cones are visible at all depths below the 
applied force. This is shown clearly by Head et. al. [30] and 
Moukarzel [31] . We use five bead-diameters because it pro- 
vides a balance between being deep enough in the pack that 
the light cones have had room to spread apart enough to be 
resolveable, and being close enough to the applied force that 
excessive averaging is not required. Even though the response 
for the single pack varied wildly throughout, on average the 
large deviations in the middle cancel, leaving stress only in 
the region where x m y. 




FIG. 4. A very small sample pack to illustrate the method of 
computing the response described in Section [IT] Three beads 
of radius R, numbered 1-3, form an equilateral triangle. They 
sit on a rough floor, so that bead 3 is touching the ground at 
a 45° angle, and bead 2 is touching the ground in two places: 
one directly below the center, and another 30° to the side. 
The six contacts are labeled A-F. Since there are six contacts 
and three beads, the system is isostatic. 



to solve this problem are available in many scientific pro- 
gramming packages. Once the contact force vector f has 
been computed, the response G can be found for each 
contact. For example, the response for contact A is lo- 
cated at (x, y) = (-R/2, —R^/3/2) relative to the applied 
force, and takes the value G = — /a(R + i?)ni 2 (ni2) a = 
( — 1/4, — V3/4), in units of the beads' diameters. The re- 
sponse at the other contacts can be computed in a similar 
manner. 
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III. SEQUENTIAL PACKS 

Previous simulations have shown the presence of light 
cones in periodic two dimensional packs built up by 
adding variably sized disks one by one onto an irregular 
floor [THl [30] • The results indicated a speed of light 
consistent with 1, meaning that the rays of force prop- 
agated out at 45° angles below a supplied point force. 
Examples of these light cones were shown in Figure [3] 
Such sequential packs are formed by taking beads one at 
a time and placing them each at the lowest point in the 
system where they will have exactly two contacts. Using 
this construction, the system has the isostatic contact 
number at each step, and the process can easily be re- 
peated until the desired system size is reached. To avoid 
crystallization, the bead radii are chosen from a polydis- 
perse distribution. Unless noted otherwise, the packs in 
this paper are bidisperse, with a radius ratio of 1.4 : 1. 
In some other similar situations, the beads have a weight, 
and gravity plays a role in the deposition process. This 
leads to the formation of local structures such as bridges 
[33) . However, here there is no gravity. The lowest point 
in the pack is determined by the starting location of the 
floor, and there are no forces on any of the beads un- 
til the point force is supplied. Only when this happens 
are there any nonzero contact forces. While in an elastic 
solid there would be a nonzero response both above and 
below the applied force, in these isostatic packs the con- 
tact forces only appear below the supplied force [3D] , and 
propagate down toward the floor in the manner indicated 
in Figure [3J 

Discovering what determines the direction of the re- 
sulting light cones was the motivating factor in this work: 
why do they go "down"? There are three obvious ways 
in which symmetry can be broken to allow for this. First, 
the direction of the applied point force provides an option 
for a preferred direction. However, this is easily shown 
not to be the deciding factor. No matter what angle 
the point force is directed in, the locations of the light 
cones do not change; the only difference is their relative 
magnitudes, as shown in Figure [5] Second, the pack- 
ing history contains a clear distinction between up and 
down: new beads are always placed above the ones al- 
ready there. However, the argument for the null stress 
law makes no reference to packing history, so changing 
this anisotropic construction should not affect the light 
cones. Third, the boundary conditions clearly provide a 
distinction between the unconstrained top and the hard 
floor. 

To separate the effects of the anisotropic boundaries 
from those of the anisotropic packing history, we con- 
struct packs in an isotropic manner, and then later break 
the symmetry by adding a floor. Because the boundary 
conditions are the same in both cases, the above argu- 
ment might lead one to expect the response to be un- 
changed. However, as detailed in Section [TV] we do not 
see light cones in this isotropic case! 
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FIG. 5. The response five bead-diameters below a point force 
in a two-dimensional sequential pack. G y , measured in units 
of bead-diameters, is the vertical response discussed in Sec- 
tion [II] and 6 is the angle of the applied force relative to the 
— y direction. This response was calculated by averaging over 
4500 realizations. Changing the direction of the applied force 
does not change the location of the peaks; only the relative 
magnitudes differ. 



IV. ISOTROPIC PACKS 

In this section we test whether an isotropic pack can 
have light cone force propagation like the sequential 
packs of the last section. Certainly no light cones are 
expected if the system gives no indication of which di- 
rection they can point. As noted above, we may define 
such a direction by adopting one feature of the sequen- 
tial pack, namely the constraining floor. We first describe 
our procedure for creating an isotropic pack in a periodic 
box. Next we show how we modify the system to create a 
floor constraint. We then discuss the effects of this modi- 
fication and argue that the isostatic character of the pack 
is not essentially altered. Finally, we discuss the peculiar 
behavior of the force response using different boundary 
modifications. 



A. Pack creation 

We follow an established method to create isostatic 
packs with no preferred direction [5]. We start with a 
periodic L x L box filled with small beads at random lo- 
cations. The bead radii are chosen according to whatever 
polydispersity we please such that the system's density 
(f) = J2i( n Ri)/L 2 is well below the critical density <p c 
where the system jams. Though (j) c varies from pack to 
pack, in two dimensional systems it is typically larger 
than 0.82 0. 

The beads are given a pairwise interaction potential 
which is harmonic if the two beads i and j are overlap- 
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ping, and zero otherwise: 
' e(l- 







' >J 
else 



< FL. 



where is the distance between beads i and j, Rij = 
Ri + Rj is the sum of their radii, and e is a constant 
that sets the energy scale. With this potential, we use a 
conjugate gradient algorithm to minimize the energy of 
the system by rearranging the bead positions. Since our 
system is fairly dilute at this point, it is possible to find 
a configuration, close to the original, where the energy is 
zero. 

Next we swell all of the beads by some percentage so 
that the packing fraction increases by a small amount 
A<j), and again perform the conjugate gradient minimiza- 
tion. We do this quasistatically so that the bead rear- 
rangements are minimal; typical A<p depend on the sys- 
tem size, but are on the order of 10 to 10~ 10 . 

Eventually the beads will be large enough that no local 
rearrangement can remove all bead overlaps, and the con- 
jugate gradient algorithm will give a system with some 
finite pressure. It has been shown empirically that at 
this transition, the system is isostatic [5]. In practice, 
due to our finite step size A^>, minor adjustments need 
to be made to get a fully isostatic pack. We first count 
the number of contacts n relative to the number of beads 
N. Note that around 5% of the beads will have no con- 
tacts. These are called floaters or rattlers, and we do 
not include them in our count of N; when the final con- 
figuration is reached, we actually remove them from the 
system altogether, since they cannot contribute to the 
force network. If n = 2N, then the system has the iso- 
static contact number and no adjustments are necessary. 
If our system is hyperstatic by a few contacts, then we 
decrease Acf> by a factor of two and reverse the sign, so 
that the beads shrink slightly instead of swelling. We 
continue to change the packing fraction until the system 
is isostatic, or it becomes hypostatic. In the latter case, 
we again take A(j> — > —A<f)/2. This is repeated until we 
either get the isostatic contact number, or A<f> gets too 
small and we give up. 

Once the system has the isostatic contact number, we 
can check to see if it is truly isostatic by counting its 
floppy modes. There are necessarily two such modes of 
T> with zero eigenvalue, corresponding to the two dif- 
ferent uniform translations available under the periodic 
boundary conditions, but these will disappear when the 
floor is added and can thus be ignored here. 

In these systems, the four directions given by the pe- 
riodic boundaries are of course the same, but this con- 
struction might allow some other direction to be differ- 
ent. Nevertheless, we refer to these systems as isotropic, 
in light of the fact that they are explicitly not anisotropic 
in the directions relevant to the sequential packs of Sec- 



tion III there is no difference between up-down and left- 
right. 

Recent work suggests that because our pack creation 
algorithm is driven purely by volume-changing processes, 




FIG. 6. On the left, a small pack created using the method 
described in Section |IV| with a dashed box indicating the 
periodic boundaries. Three beads have been labeled for com- 
parison with the image on the right, which denotes the same 
pack with a floor. Beads which cross the horizontal boundary 
(e.g. A and B) are treated as being on top of the pack. This 
means that the force balance equations for these beads do not 
include the contact with bead C. However, the force balance 
equation for C remains unchanged; what was formerly the 
contact from C to A is now a contact from C to the floor. In 
this way, both the total number of beads in the pack and the 
total number of contact forces remains unchanged. 



it will be stable under uniform compression but not under 
volume-preserving shear forces. This can result in sys- 
tems with a negative shear modulus [34 . Furthermore, 
due to the extra degree of freedom allowed by boundary 
distortions, the number of contacts required for a finite 
system to jam may not be the same as the isostatic num- 
ber discussed earlier. In this case, rather than using the 
lack of floppy modes as the determining characteristic for 
jamming, the existence of nonzero bulk and shear moduli 
can be used instead [35]. However, all of the systems that 
we study are modified versions of these periodic packs, 
where at least one of the boundaries will be rigidly fixed 
in place, removing its degree of freedom. We find that 
counting floppy modes is sufficient to make our pack sta- 
ble under the forces we apply. 

To break the symmetry, we now create a bumpy floor 
to resemble the sequential case. This is done by tak- 
ing some arbitrary horizontal line across the pack, and 
designating all beads it crosses to be on the boundary. 
Each boundary bead serves both as a floor for the beads 
above it to rest on, and also as a bead sitting on the top 
of the pack. This is accomplished by keeping the force 
balance equations the same as they would be otherwise, 
except that for beads on the boundaries, contacts with 
non-boundary beads above them are not included. This 
is explained in more detail in Figure [6j 

This modified system has the same number of beads 
and contacts as the old one, and thus still globally satis- 
fies the isostatic contact number. However, there are now 
fewer force balance equations on the bottom, since none 
of them are used for the floor. Because the number of 
contacts did not decrease, there is a locally hyperstatic 
region of the pack here. On the top, we have a corre- 
sponding hypostatic region. This causes two problems. 
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FIG. 7. On the left, an image of one of the hyperstatic null 
vectors of the contact matrix M. M becomes singular when a 
hard floor is added using the method shown in Figure [6] The 
lines between beads represent contact forces, with color de- 
noting sign and thickness denoting magnitude. Even if some 
multiple of this nonzero pattern of contact forces is put on 
the beads, there will still be zero net force on all beads. On 
the right, we see that the null vectors are localized near the 
boundary. For several different pack sizes, we look at the 
null vectors formed when the hard floor is added. Each null 
vector is scaled so that its L2 norm is 1, and then the aver- 
age magnitude of all contact forces at a particular height is 
computed. The average is taken over all null vectors created 
in each of several different floor configurations. The differ- 
ent floor configurations were chosen by putting the horizontal 
boundaries that determine them at several heights, separated 
by a distance of two bead-diameters. No matter the size of 
the pack, the strength of the null space decayed following an 
exponential whose decay length was between two and three 
bead-diameters. 



First, the existence of a hyperstatic region means that 
the contact matrix M has at least one null vector involv- 
ing the contacts on the bottom of the pack. This is a 
nonzero pattern of contact forces which sum to zero net 
force on every bead. Thus there is no unique solution 
to our equation Mf = F, as arbitrary multiples of these 
null vectors can be added to f to give the same F. This 
problem can be solved by finding the f which minimizes 
the compressive energy. As discussed in Section [TTJ this 
is done by choosing the f for which |f| 2 is a minimum. 

Though we give a prescription for choosing a particular 
f, we do not actually expect this choice to matter in the 
bulk of the pack where our measurements are made. Fig- 
ure [7] shows that the magnitudes of the contact forces in 
the hyperstatic null vectors decrease exponentially with 
distance from the boundary, with a decay exponent on 
the order of a few bead-diameters. Thus while multi- 
ples of these null vectors can be added to any solution, 
only the contact forces near the floor are affected. More- 
over, the number of null vectors created by the process 
of adding a floor scales with yN, the linear dimension 
of the system. As the number of beads N in the system 
increases, these null vectors become a vanishingly small 
fraction of the number of modes in the system. This 
scaling is shown in Figure [8j 

The second problem caused by the creation of a null 
space of M is more worrisome. Because the system is 
hypostatic at the top, M is rank-deficient. This means 
that the space of all vectors Mf has a smaller dimension 
than the space of possible body forces F. Thus there is 



FIG. 8. When a hard floor is added to the system using the 
method shown in Figure [6] the contact matrix M becomes 
singular. In (a), we show that the average number of null 
vectors of M scales with the linear size of the system; that 
is, with the square root of the number of beads N. Thus as 
the system size is increased, a vanishingly small percentage 
of the modes of M will have a zero eigenvalue. In (b), we 
show the probability that a certain null space dimension was 
measured, given the size of the system. The lines show a 
Poisson distribution fitted using a mean equal to the average 
null space dimension for that system size. 

only an infinitesimal probability that some F will actu- 
ally have a contact force network f which can support it; 
the system is unstable. In particular, it is in general im- 
possible to find any solution to the problem of supplying 
a point force. 



B. Applying forces to a local cluster of beads 

One way to guarantee that the system will be stable 
under an applied force is to modify the force so that it is 
in the image of M. This can be done by applying external 
forces to several beads in a small localized cluster, rather 
than to a single bead. If the null space of M has dimen- 
sion k, then the vectors Mi span a space of dimension 
2N — k with N the number of beads in the system. Thus, 
some combination of k linearly independent forces can be 
added together to give an F that is in this 2N — k dimen- 
sional image of M. We can take the linearly independent 
forces to be, for example, vertical and horizontal forces 
on k/2 beads in a small region. Furthermore, if we use 
k + 2 independent forces, then we can specify the mag- 
nitude and direction of F as well. This lets us supply 
a force similar to the point force used in the sequential 
packs, though it is spread out into a slightly larger area, 
as shown in Figure |9| Since k scales like y/N, this cluster 
becomes a better approximation of a point force as the 
system size is increased. Strangely, it is typically pos- 
sible to obtain an F in the image of M with far fewer 
than k + 2 independent forces. This must be due to some 
special feature of M's structure, as this is not true for an 
arbitrary matrix. We have unfortunately been unable to 
determine what this feature is. Nevertheless, using this 
method, the contact forces can be found for these packs, 
allowing the computation of a well-defined G. 

As shown in Figure [3^,, for an individual pack the re- 
sponse looks quite wild, and many realizations need to 
be averaged together for the continuum behavior to be 
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FIG. 9. Instead of supplying a force to a single bead, we 
can apply forces as shown to a small cluster of beads. In 
this case, the pack's contact matrix had a 16 dimensional 
null space, so we applied x and y forces to 8 beads in the 
circle indicated in on the left, to give the net downward force 
indicated. On the right, we show the individual forces applied 
to the beads. Their magnitudes were chosen so as to minimize 
|f| 2 , as discussed in Section |ll| 



apparent. To get an idea of how many realizations are 
needed, we can construct a function c, defined by 
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FIG. 10. The function c of Equation [IT] for both sequential 
and isotropic packs, a distance of five bead-diameters below 
the location of the applied force, c depends slightly on the par- 
ticular choice of realizations used to compute G^' and ■ 
The error bars show the spread of values obtained over 25 
different choices. For sequential packs, c decreases as 1/y/R. 
For isotropic packs, it shows no sign of decreasing, even over 
three decades. 



fAT , Y, Xl \G[ NR \x l ,y,)~G[ N -\x l ^)\ 
c(N R ,y ) = t^—, (17) 



C. Modifying the boundaries 



where g[ Nr " > and G^^ are each responses averaged 
over different sets of Nr realizations. For any particular 
?/o, this indicates how closely the two different functions 
match. If c(Nr,uo) is small, then both sets of Nr re- 
alizations average to something similar, and Nr is large 
enough to see the continuum behavior. 



For the sequential packs discussed in Section III 



scales like 1/^/Nr, as could be expected. However, as 
seen in Figure [TUJ c does not decrease with Nr for the 
isotropic packs, indicating that averaging does not give 
any convergence; there seems to be no continuum be- 
havior that would allow one to speak of a coarse-grained 
stress field. 

This lack of averaging is possible due to the distribu- 
tion of contact force sizes. Disorder in isostatic systems 
causes multiplicative rather than additive noise. Sys- 
tems with this sort of noise do not average well: the 
Green's functions have magnitudes that vary in a range 
that grows exponentially with distance from the applied 
force, and follow a power law distribution P(G) ~ G~ a , 
with a close to one [31]. Thus in any collection of packs, 
the magnitudes of the contact forces can vary exponen- 
tially, and a small number of packs with large forces end 
up dominating the average response [36 . At any fixed 
depth below the forcing point, Newton's Laws assure us 
that the total force summed across the entire pack width 
is equal to the supplied force. However, the higher mo- 
ments diverge, so any averages over a small portion of the 
pack, rather than the whole width, can fluctuate wildly, 
making it impossible to measure a continuum response. 



The null modes of the previous section appear in the 
process of adding the floor. Instead of using a straight 
horizontal line to determine the floor beads, more care 
can be used to create a floor which does not lead to any 
null modes. To find such a floor, we use a simulated 
annealing algorithm to minimize the dimension of the 
null space by making small deformations to a starting 
boundary, as demonstrated in Figure |11| 

The simulated annealing algorithm starts by creating 
a directed graph Q whose nodes are the beads in the pack 
and whose edges are the contacts, directed from the left 
to the right. The periodic boundary conditions in the 
horizontal direction allow any path on this graph from 
a bead back to itself to be treated as a floor. Making 
the graph edges point from right to left will also work as 
long as the choice of direction is consistent; the idea is to 
create a floor that spans the entire horizontal length of 
the system instead of immediately looping back on itself. 

Given any floor, the size of the null space of M can 
be computed. This can be used as some sort of effective 
energy £ that we are trying to minimize. Neighboring 
floors are chosen at random by picking two nearby beads 
i and j in the current floor and finding the shortest path 
between them on the graph Q' , which is Q with all edges 
on the current path between i and j removed. On rare 
occasions the new path will be much longer than the 
original. This is typically a sign that the path wraps 
around the periodic boundaries multiple times, and so 
such paths are discarded. The new path will have some 
energy £' associated with it. We keep this path with 
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FIG. 11. On top, a floor made from a straight horizontal 
line across a pack. The null space for this system is six di- 
mensional. A combination of the null vectors is shown using 
black and gray lines. Each line represents a contact force be- 
tween the two beads, whose magnitude is proportional to the 
line width and whose sign is indicated by the color. The null 
vector shows a way to have nonzero contact forces, but nev- 
ertheless have zero net force on every bead. In the middle 
image, the center of the floor has been deformed to eliminate 
part of the null space. The resulting system's null space is 
only four dimensional. On the bottom, the floor has been 
further modified so that there is no null space; the contact 
matrix is now invertible. 



FIG. 12. The density of states for a 2000 bead pack before 
and after a hard ceiling has been added using the method de- 
scribed in Section|V] The general shape is unaffected, and the 
low frequency plateau expected in an isostatic pack is evident 
even after the change. The only significant difference is in 
the first data point, which counts the number of modes be- 
low some small Aui. This difference is due to the null modes 
discussed in Section |IV| which vanish when the hard ceiling 
is added. As noted in the text, the fraction of modes af- 
fected this way becomes vanishingly small as the system size 
increases. 

One cannot speak of any continuum stress field, let alone 
one which obeys a null stress condition. 



probability 



V. STABILIZING THE ISOTROPIC PACKS 



P 



-(s'-syr i{£ i >£ 
otherwise 



(18) 



The effective temperature is chosen to go like T ~ 
l/log(s), with s the number of steps taken so far. This 
temperature function can probably be improved, but the 
process finds states where £ = fairly rapidly as it is. 

When a floor with £ = is found, the algorithm will 
remember it and continue to run until some set number of 
steps, typically on the order of 10 3 , is reached. It is rea- 
sonable to wonder whether the existence of floors which 
create no null space is dependent on having a finite-sized 
pack. This turns out not be the case, however. In fact, 
the larger the pack, the more likely the existence of such a 
floor. This is shown in Appendix[Bj Thus, using a variety 
of starting points, many floors that do not cause a null 
space can typically be found. As discussed in Section \U\ 
we expect the lowest modes of T> to be the most im- 
portant for determining static force distributions. From 
the set of possible floors, we choose the one whose low- 
est eigenvalues have magnitudes similar to those found 
in sequential packs of the same size. 

Even in this case, where the important eigenvalues of 
the isotropic packs are similar to those of the sequential 
packs, the response varies too wildly to give an average. 



In the previous section, we showed that periodic 
isotropic packs can be modified to have a floor. Though 
this creates some issues by making M singular, these 
problems can be overcome to give a well-defined force 
response for a given pack. However, the resulting forces 
vary so wildly from pack to pack that it is impossible to 
infer a continuum average. Here, we make a more drastic 
change to the pack in an attempt to stabilize it and cre- 
ate a continuum response. We do this by adding a hard 
ceiling to the pack, in addition to the hard floor. Instead 
of imposing force balance constraints for the boundary 
beads (such as beads A and B in Figure [6| we simply 
allow them to absorb any force. 

This system is no longer isostatic, as it has excess con- 
tacts near both the floor and the ceiling. However, the 
bulk of the pack is unaffected, and remains isostatic. 
Adding the ceiling does not greatly affect the density 
of states; aside from removing the zero eigenvalues, the 
overall shape of the plateau is not altered. This can be 
seen in Figure 12 As noted in the previous section, the 



number of zero eigenvalues of the pack without the hard 
ceiling goes like v^V, so the number of affected modes is 
an asymptotically small fraction of the total number. 

Unlike the case with a free top, the responses now av- 
erage nicely. However, instead of forming light cones, the 
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FIG. 13. A comparison of the response five bead-diameters 
below the forcing point for 5584 isostatic packs with the same 
number of hyperstatic [Z ~ 5.6) packs, when hard ceilings 
are added. The hyperstatic packs were formed by adding 
contacts to the isostatic ones between all pairs of nearest- 
neighbor beads. The responses are in arbitrary units, normal- 
ized to the maximum value. Both pack types have the same 
functional form. In both cases, the response more closely re- 
sembles an elastic solid, with a single peak below the forcing 
point, than it does the sequential packs. 



stress now peaks directly below the applied force. This 
is shown in Figure |13| The response of our isostatic pack 
with slightly hyperstatic boundaries now resembles that 
of a system which is hyperstatic everywhere. That is, 
the average force response behaves more like it would in 
an elastic solid [TH [37] . The isotropic packs must then 
in some sense be stronger than the sequential ones, since 
one can restore elastic-like behavior to them with much 
less modification. 

It may be that this new behavior in the bulk is due 
to the hyperstaticity at the boundaries. However, this 
should not be the case. Even when the average number 
of contacts per bead in some pack exceeds the isostatic 
number by some amount Az, it should only behave elas- 
tically on length scales larger than t* ~ 1/ Az [TTJ [T^] . 
There is evidence that an ensemble average of many such 
packs with a small Az will nevertheless show elastic-like 
behavior on smaller length scales when the contacts are 
added randomly throughout the system [5]. However, 
this may depend on the ensemble average being able to 
distribute the resulting elastic energy stored in the extra 
bonds uniformly throughout the pack. In our system, the 
excess contacts are all near the boundaries. This means 
that not only does the fraction of extra contacts Az —> 
for large packs, but also that any sort of elastic behav- 
ior should be limited to the contacts there. Unless the 
coarse-graining is done on length scales larger than the 
distance between the ceiling and the floor, there is no 
way for the extra contacts to affect the behavior in the 
bulk. We can test this by giving sequential packs the 
same floors and ceilings as these isotropic ones, which 
yields similar departures from isostaticity. This does not 




FIG. 14. On the left, a plot of the stress a yy in isotropic packs 
with hard ceilings, obtained by averaging over 15254 realiza- 
tions. For comparison, on the right is a yy for the sequential 
packs when hard ceilings have been added. In both cases, 
the forcing point is in the center, and is directed down. The 
light cones are still clearly visible for the sequential packs, 
indicating that the presence of the ceiling does not cause the 
elastic-like behavior seen in the isotropic packs. 



create a stress profile with a single peak. Instead, the 
light cones are visible going to both the floor and the 

This further indicates 



ceiling, as shown in Figure 14 



that it is the isostatic bulk of the isotropic packs which 
cause the elastic-like response, and not the boundaries. 

The lack of light cones in the isotropic packs with hard 
ceilings is an indication that the null stress law expressed 
in Equation [2] is not satisfied. Previous verifications of 
the null stress law have been done by applying a uniform 
load across the top of the entire pack, and then measur- 
ing the average &}■(.■ By varying the direction of the load, 
the <jm can take on a range of values. Head et. al. [30 
found that a xx /a yy was essentially constant as long as 
Cxy/cyy 0.3. This was consistent with the null stress 
law of Equation[2]with /i = 0, though it broke down when 
the shear became large enough that the left-right symme- 
try in their systems was broken. However, this approach 
is not sufficient to distinguish null stress behavior from 
elastic-like behavior. 

To see this, say that the load applied evenly across the 
top of the Lx L pack is F = (F x ,F y ). Because there can 
be no x dependence in this case, Equation [T] reduces to 
9 y cr xy = d y a yy = 0, implying that those components of a 
also have no y dependence. Their values can be found at 
the top: a xy — F x /L and a yy = F y /L. If the null stress 
condition holds with fi = 0, then as F varies, cr xy /a yy 
changes, but cr xx /a yy = r/ remains constant. This is in- 
distinguishable from what happens in a two-dimensional 
elastic solid. There, the continuity equations give the 
same values for a xy and a yy . To find a xx , note that the 
strain satisfies 



^xx — dx^x 



K(a xx — a yy ) + G{a xx + a yy ) 



(19) 



AGK 

where K and G are the bulk and shear moduli, respec- 
tively. The bulk modulus measures the system's resis- 
tance to uniform compression, and the shear modulus is 
the ratio of shear stress to shear strain at constant vol- 
ume |15j . Since there is no x dependence, Equation 19 
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tells us that u xx /a yy 
is the Poisson ratio. 



= is, where v = (K - G)/(K + G) 
Again, as F varies, cr xy /a yy will 



change, but a xx ja yy will remain constant. 

To distinguish between the two situations, we do not 
supply a load across the entire width of the pack. In- 
stead, we supply a point force as elsewhere in this paper, 
and then use Equations [12] and [13] with w = 4 bead- 
diameters to find the average stresses at several random 
locations. For each such area, we then compare <J xx j a yy 
with a xy /a yy . If the null stress condition holds, then 
&xx/vyy should remain essentially constant, no matter 



the value of a xy . Figure 15 shows the results for both 
isotropic packs with hard ceilings and sequential packs. 
For comparison, we also show the stresses from a peri- 
odic two-dimensional elastic solid with hard floors and 
ceilings. The elastic stresses were computed using the 
results of Leonforte, et. al. [38], who extended the 
work of Serero et. al. [55]. The relevant formulas are 
quoted, using our notation, in Appendix[C] We find that 
in the isotropic packs with added hard ceilings, there is 
clearly no linear relationship between the components of 
the stress tensor. The sequential packs, on the other 
hand, seem to be more in line with the null-stress behav- 
ior we anticipated. However, this test is only intended 
to demonstrate the clear contrast between the sequential 
packs and the isotropic ones; we have not unambiguously 
identified the source of this contrast. 



A. Effective moduli of isotropic packs 

Since the null stress relationship between the compo- 
nents of the stress tensor is not satisfied for the isostatic 
packs with hard ceilings, we can ask what other consti- 
tutive equation might describe them. As noted above, 
the response in many ways resembles what one would see 
in an elastic solid. Thus, we can compare the results to 
linear elasticity and attempt to find effective moduli. 

To compute the moduli, we supply a vertical force 
F = Fy to each of the beads that are intersected by a 
horizontal line through the middle of the pack. This ap- 
plied pressure will cause a discontinuity in the yy stress. 
Above the line it will take some value a yy , and below the 
line it will be —cr vy . In a two dimensional elastic solid, 



the stress-strain relationship is 



2G 



-u 



(20) 



where repeated indices are summed jTH]. The factor of 
1/2 comes from working in two dimensions rather than 
three. We do not apply any shear with our force, so <J xy , 
and thus u xy , are both zero. Furthermore, the applied 
force does not allow any x dependence, so u xx — d x u x — 
0. This leaves us with 



Ku y 
Ku v 



Gu y 
Gu y 



(21) 
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FIG. 15. The relationship between the different components 
of the stress tensor, computed using Equation |12[ with a 
coarse graining function of width w = 4d. Each data point 
gives the ratio of the stresses at a random location. The 
black diamonds were obtained by averaging over 12474 real- 
izations of the sequential packs of Section [111] and the squares 
show the results for isotropic packs with added hard ceil- 
ings, as described in Section [V] All packs were made using 
N — 500 beads. For comparison, the crosses show the re- 
sults for an elastic medium, computed using the formulas in 
Appendix [C] The measurements fall in a narrow annulus be- 
cause the computed values are the average over an area. The 
elastic medium was given a Poisson ratio of v = 0.9, which 
is the value estimated for the isotropic packs of this size (see 
Figure [l6]) using the methods described in Section [V] The se- 
quential packs approximately satisfy the null stress condition 
of Equation[2] though a- xx /a yy evidently varies by up to 25%; 
the larger values come from points nearer to the source. The 
isotropic packs, on the other hand, clearly do not. Though 
there are some quantitative differences with the elastic media, 
the general form is the same. While not evident from this fig- 
ure, the isotropic packs have two distinct values of a xx ja yy 
where <J xy /tj yy is zero, just as in the elastic case: one near 
<Txx/vyy ~ 0.4, and another around <J xx /a yy ~ 10. This lat- 
ter ratio occurs at points far from the applied force, where 
o yy is close to zero. Small numerical fluctuations in a yy thus 
cause large deviations in the computed ratios, so such points 
are not plotted in the figure. 



Dividing gives 



@xx 

a yy 



K-G 
K + G 



(22) 



Thus <J xx /<J yy is the Poisson ratio, just as it was in the 
case of the unconstrained top discussed earlier. 

To determine the bulk and shear moduli separately, 
we measure the energy stored in the pack. The energy 
density in a two-dimensional elastic solid can be written 
in terms of the strains as |15j 



-Ku 2 ^ + G[u a p 



(23) 



with the factor of 1/2 again coming from the dimension- 



ality of the system. From above, u x 



and 
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= <j yy /(K + G). The energy stored in the entire 



pack is therefore 



£ = L 2 e = 



L 2 



uu 



2 K + G 



(24) 
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to compute 



We can also use £ = ^f T f from Equation 
the energy stored in response to our force! - Equating the 
two gives 



K + G = k 1 ^ 



(25) 



Combining this with Equation [22] lets us solve individu- 
ally for the moduli: 



K= kL2c vv( 
2 |f| 2 V 

= k L 2 a 2 yy 

2 Ifl 2 



1 



a 



XX 



(7 



1 - 



yy 

^ XX 



(26) 



a 



'!'! 



Note that as K gets large compared to G, the formula for 
G becomes less accurate. An alternative form of forcing, 
F = Fx, would provide a shear which is not dominated 
by the effects of the bulk modulus. 

Since a xx and a yy are constant on either side of the 
applied force, they can be computed using Equation |T2"| 
with a <& that encompasses either the top or bottom half 
of the pack. We expect the area away from the bound- 
aries and the applied force to be better-behaved, so we 
can imagine an L x 3L/10 area A which includes the 
half of the pack between the applied force and the ceil- 
ing, excluding anything within L/10 of cither of those 
boundaries. Then 
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A perfectly ordered isostatic lattice can have finite 
moduli, but in disordered isostatic systems, they ap- 
proach zero as the system size increases. Moukarzel [ID] 
has presented numerical data suggesting that the form 
of the moduli depends on the packing symmetries. For 
isotropic isostatic networks with periodic boundary con- 
ditions in all directions, he found that they decrease like 
AN~ B , where B depends on the pack construction and 
disorder. However, for packs constructed from a hard 
floor with a preferred direction, the decay is exponential. 
He argues that this is caused by the multiplicative noise 
discussed in Section ftVBI 

The systems we consider here share the isotropic bulk 
of the power-law systems, and the preferred boundary di- 
rections of the exponential systems. Figure [16] shows that 
moduli follow the power-law behavior of the isotropic 
packs, with exponents of Bk = 0.21 ± 0.01 for the bulk 
modulus and Bq = 0.85 ± 0.07 for the shear modulus. 
With the added hard ceiling, our system does not dis- 
play the exponentially large contact forces, supporting 
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FIG. 16. A log-log plot of the bulk and shear moduli of the 
isotropic packs with hard ceilings and floors of Section W\ as a 
function of the number of beads N which make up the pack. 
The moduli are normalized by the spring constant k of the 
harmonic force between contacting beads. Each data point is 
the average of at least 100 different realizations of the same 
size. The dashed line is a fit to the form AN~ B , with Bk = 
0.21 ± 0.1 for the bulk modulus, and B G = 0.85 ± 0.07 for the 
shear modulus. The beads in the packs have a bidispersity of 
1.4 : 1. Moderate changes to the bidispersity do not affect 
the scaling. 



the claim that it is the multiplicative noise and not the 
boundary conditions which cause an exponential decay. 
An alternative explanation given for this elastic collapse 
in Moukarzel's systems with average contact number of 
exactly Z = 4 is that such systems are not actually 
jammed, due to an extra degree of freedom which comes 
from the system size |35j . Since the systems are not 
jammed, vanishing moduli are expected. However, our 
systems are over-constrained by more than one contact, 
so the single degree of freedom used in this argument 
does not seem to explain our decay, though it is possible 
that we have simply not gone to large enough N . 

Because the shear modulus decays faster than the bulk 
modulus, the Poisson ratio approaches 1 in the limit of 
infinite system size. This is the expected value of v: in 
disordered packs with small excess contact number Az, 
the bulk modulus is proportional to the effective spring 
constant k, whereas the shear modulus is proportional to 
kAz [5]. In two dimensions then, v = (K—G)/ (K+G) —> 
1 as Az — > 0. 

These systems with added hard ceilings show that it 
is possible to have a continuum response in an isostatic 
system with no light cones. The stress can instead fol- 
low a pattern that can be described using linear elas- 
ticity, with bulk and shear moduli which exhibit power 
law decays when the system size increases. The presence 
of elastic-like behavior indicates that light cones do not 
arise from simple isostaticity. More specific geometric 
concerns must also play a role. These are discussed in 
the next section. 
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FIG. 17. On the left, we show the contact angles for a bead. 
They are defined to be the angle between each contact and a 
line pointing up from the floor, as shown by the dashed lines. 
Note that if for some bead we measure an angle 9, then for 
the contacting bead, we will measure tv — 9. The distribu- 
tion is therefore symmetric about n/2. On the right, we see 
the probability density function P(8) found by computing all 
9 values in 100 packs, for both the sequential and isotropic 
cases. While the isotropic packs have no preferred angle, the 
sequential packs have strong peaks. Due to a preference for 
horizontal rather than vertical contacts, the average contact 
angle is slightly to the side of these peaks, close to (9) ~ 45°. 



VI. THE CONTACT ANGLE DISTRIBUTION 
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FIG. 18. The probability density function (PDF) of contact 
angles P{9) for a pack with a large bidispersity. The radii of 
the large beads are 3.5 times bigger than those of the small 
ones. The solid line shows the PDF using all of the contacts 
in 100 packs, and is symmetric about 7r/2. The two other 
curves show the PDF for only the contacts of beads of a par- 
ticular size, weighted by the fraction of the total number of 
contacts in the pack with beads of that size. These two curves 
sum to the PDF for the entire pack, but are not individually 
symmetric. 



Though many properties of sequential and isotropic 
packs are similar, their contact angle distributions are 
quite different. For a pair of beads, we can look at the 
angle 9 which the line between their centers forms with 
respect to a vertical line from the floor. By doing this 
for every pair of contacting beads in many packs, we get 
a probability distribution P{&). 

In our isotropic packs, P(ff) is a constant. However, 
in the sequential packs, there are peaks near 45° and 
135°, as shown in Figure [P7) This is the same angle at 
which the light cones propagate. The correlation of light 
cone direction with bond angle suggests a more prosaic 
explanation for the light cones than the formal null-stress 
picture. Light cones may arise in sequential packs simply 
because each bead is sitting on top of two others, which 
on average will be in the directions the light cones go. 

We can attempt to test this by varying the distribution 
of the bead sizes in sequential packs. This will modify 
the probability density function P{9) of the contact an- 
gle distribution, which can be compared to the positions 
of the light cones. In addition to a few polydisperse sys- 
tems where the radius of each bead is chosen from the 
flat distribution in the range [R, aR], we also use several 
bidisperse systems where the radius is either R or aR. 
In the bidisperse cases, half of the beads are given each 
radius, so that the larger beads occupy more total area. 

For the more extreme bidispersities, P(9) no longer 
has the form of two clean peaks, as shown in Figure [18] 
However, most of the deviation from that form is due to 
the small beads, which end up not being important for 
the force transmission. Many of them serve only to fill 
in the spaces between large beads, and do not experience 
any contact force; while typically about 80% of the large 



beads in a pack experience some nonzero contact force, 
only about 15% of the small beads do. Thus when com- 
puting the average contact angle for comparison with the 
light cone direction, we use the PDF for the downward 
propagation from the large beads only, since they are the 
ones dominating the contact force network. 

Table [I] compares the average 9 in the range [0, 7r/2] 
with the angle tp the light cones make with the vertical, 
for a range of radius distributions. While 9 is consistent 
with tp for the cases shown, the range of possible angle 
distributions that can be obtained from sequential packs 
in this way is quite limited, and any shifts in the peak 
positions are within uncertainties. 

To generate more extreme changes in the contact angle 
distribution, we look at packs formed from isostatic lat- 
tices. For each of the base lattice configurations shown in 



Figure 19 we perturb the position of each bead by mov- 
ing both its x and y coordinates a random amount cho- 
sen from a normal distribution of width f3R. This causes 
some beads to overlap, but we ignore this and use only 
the contact angles to construct M and find the response. 
This effectively treats the system as point particles con- 
nected by springs, and the positional randomness serves 
merely to modify the rest lengths and angles of these 
springs. 

Floors are added to these packs as shown in Figure[6]for 
the isotropic case. The diagonal lattice is not so interest- 
ing; the average contact angle is still 45° , and this is also 
the angle of the light cones. However, the kagome lat- 
tices do show different distributions, which are depicted 
in Figure [20} In the orientation of Figure [T9p, there are 
three peaks in the contact angle distribution. One corre- 
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Type 


a 


(9) 


V? 




1.4 


43.4 ± 0.6° 


40 ±3° 




2.0 


41.7 ± 0.6° 


38 ±4° 


Bidisperse 


2.5 


39.4 ± 1.0° 


40 ±5° 




3.0 


42.8 ± 1.1° 


40 ±5° 




3.5 


45.8 ± 1.1° 


38 ±4° 




1.1 


44.0 ± 0.4° 


44 ± 2° 


Polydisperse 


2.0 


43.0 ± 0.4° 


40 ±3° 




3.0 


42.5 ± 0.6° 


40 ±5° 




6.0 


43.5 ±0.6° 


39 ±4° 



TABLE I. A comparison of the light cone positions tp with 
the average contact contact angle (8), for several different 
pack types. All of the packs used here were created using the 
sequential method of Section |III| For the bidisperse packs, 
half of the beads were given radius R, and the other half had 
radius aR. For the polydisperse packs, the radii were chosen 
randomly from the uniform distribution [R, aR] . In the bidis- 
perse packs, the disorder in the pack initially increases with 
a, creating wider and less defined light cones. As a increases 
further, the small beads become insignificant compared to the 
larger beads, and the system becomes more ordered; the small 
beads play only a minimal role in both the overall structure 
of the pack and the force transmission. This is discussed in 
Figure |18| The polydisperse peaks tend to be more clearly 
defined even at large polydispersities, as there is no sharp 
separation of beads based on size. 









FIG. 19. A diagonal lattice (a) and two different orienta- 
tions of a kagome lattice (b and c). All of these arrangements 
have the isostatic contact number. When they are used, ran- 
domness is added in the form of modifications to their beads' 
positions. Bead displacements in the x and y directions are 
selected at random from a normal distribution of width PR 
for R the bead radius and ft some positive number. This 
randomness serves not only to spread out the contact angle 
distribution, but also to remove the null modes inherent in 
the unperturbed kagome lattices. 





sponds to horizontal contacts, and the others are around 
(9) = 30.0 ± 1.4°. The light cones in these kagome packs 
are located at <p = 30 ± 3° . In the other orientation, one 
of the peaks in the contact angle distribution corresponds 
to vertical contacts. The origin of this is obvious from 
Figure [T9p. For each realization used in the average, the 
applied point force was exerted on a random bead. Thus 
a large proportion of the forcing points were located on 
one of these vertical struts, since that is where most of 
the beads are. In this case the simplest direction for the 
force to go is straight down the strut. This gives the 
large peak in the response corresponding to light cones 
at ip = 0. The other light cones are at <p = 60 ± 3°. 
These correspond to the average value of the other peaks 
in P(d), at {6) = 59.8 ± 2.9°. 

Thus in the isostatic lattices we also find that the light 
cone angle follows the average contact angle. Given the 
ordered nature of lattices, this is not as surprising as it 
was in the much more disordered sequential packs. Ex- 
perimental studies have shown that the forces in regular 
granular lattices such as hep and fee crystals indeed move 
in rays following the preferred lattice directions, without 
any need for isostaticity |41| . When a small amount of 
disorder was added to such systems, there was a tran- 
sition to a broad central peak after a few layers. Here, 
we see the light cone-like behavior at all depths in our 
isostatic lattices, even with significant disorder. Admit- 
tedly, our packs are small enough that we can only verify 
this to a depth of about 20 bead diameters. However, 
since we also see the same light-cone like behavior in the 
disordered sequential packs which have no obvious under- 
lying lattice structure, we believe that in this case, it is 



FIG. 20. A comparison of the contact angle distribution with 
the light cone direction for the two different orientations of 
the kagome lattice. In both cases the positions of all beads 
have been perturbed by an amount chosen from a normal 
distribution with a width of one fifth of a bead radius, (a) 



shows the results for the orientation depicted in Figure 19 b). 
The two non-horizontal peaks in P(0) wee at (8) ~ ±30°, 
and the light cones clearly propagate at the same ang le, (b) 



shows the results for the orientation depicted in Figure 19 c). 
There is a peak in P(8) at 0/n, which gives the strong peak in 
the center. The other peaks in P{8) are at (8) ~ 60°, which 
agrees with the angles of the smaller light cones. The middle 
peak is larger because forcing points were chosen at random, 
and most of the possible choices were on one of the vertical 
struts. Both sets of response data were obtained by averaging 
over around 450 realizations of N = 2340 beads. 



the isostatic properties rather than the lattice ones which 
cause the light cones. Nevertheless, it would be nice to 
have an example of a system with both a substantially 
different (9) from what is seen in sequential packs, and 
no clear lattice structure. 

We attempt to find such systems by taking the disor- 
dered sequential packs and modifying their contact net- 
works to give some other target P(9). We do this by 
constructing a list of the nearest neighbors of each bead. 
We then try to randomly remove one existing contact 
from the pack, and then add another one from the list 
of nearest neighbor pairs that are not in contact. If this 
move makes the contact matrix M singular or takes our 
P{9) further from the target distribution, then we reject 
it and try again. Otherwise, we keep that move and make 
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FIG. 21. The original contact angle distribution for a single 
sequential pack, along with the target (flat) distribution that 
we want to achieve by changing contacts according to the 
algorithm described in Section|Vl] The final distribution after 
the angle modification algorithm has completed still shows 
traces of a peak, but is on a whole much flatter. 



another one. This can be continually repeated. In prac- 
tice, we continue attempting these contact swapping op- 
erations until less than one percent of them are accepted. 
This lets us get very close to the target distribution, as 



shown in Figure 21 



If we start with a sequential pack and its requisite 
peaked P(9) and then change it to have a flat distribu- 
tion, the light cones disappear; as in the isotropic packs, 
the force distribution becomes too irregular to charac- 
terize. However, if we start with an isotropic pack and 
make the angle distribution peaked, we do not create 
light cones. Likewise, if we take the sequential pack with 
the flattened distribution from before and then try to 
restore the original contact distribution, the light cones 
do not reappear. There is evidently some other aspect 
of the pack which is being destroyed by this procedure. 
This other aspect must be due to the isostatic nature 
of the pack, since in a hyperstatic pack with a similarly 
peaked angle distribution, the response is a single large 
peak directly under the forcing point, rather than light 
cones. It seems that both isostaticity and a peaked P{9) 
are necessary, but insufficient, to cause light cones. 



VII. CONCLUSION 

In summary, we have shown that the set of jammed iso- 
static systems is richer than formerly thought. There ap- 
pear to be at least three different types of force propaga- 
tion possible: the previously observed light cones which 
appear in sequential packs; a more elastic-like response 
in stabilized isotropic systems; and a non-continuum be- 
havior where no average stress can be found. The light 
cones are not visible in all isostatic packs because the null 
stress condition holds for only a subset of them. The pre- 
vious observations of the null stress law seem to depend 



not only on the isostatic nature of the packs, but also on 
some further feature which appears related to anisotropy 
in the contact angle distribution. 

Our simulations were motivated by the desire to probe 
the marginally jammed state in ways which are diffi- 
cult for experiments. Many experiments which consider 
the response to point probes use macroscopic particles 
such as photoelastic disks or grains of sand. Such sys- 
tems tend to be hyperstatic rather than isostatic, have 
frictional forces, and do not allow for tensile contacts 
[2"31 1571 H2l 22]. Other experiments look at crystalline, 
rather than disordered, packings [UJ. Still others are 
more interested in the resulting particle rearrangements 
[H] or the critical scaling close to isostaticity [35] than 
they are in the force response. 

The existence of friction in experimental systems is 
one barrier to performing real-world tests of our results. 
In the presence of these extra forces, there is no longer 
a contact number directly correlated with the onset of 
pack stability [35] . This may affect the force propagation 
at isostaticity. Frictional forces may also decrease the 
range of length scales over which hyperbolic behavior can 
be observed in packs which are only slightly hyperstatic. 
Their presence tends to make the response as a whole 
more elastic-like |47j . This barrier can be overcome by, 
for example, using a colloidal system where the beads 
do not directly touch, eliminating any contact friction 

mama. 

We allowed tensile contacts in our systems based on 
the discussion in Section [TT] which noted that the pres- 
ence of light cones in sequential packs is not affected by 
restrictions on the sign of the contact forces. However, 
the sign may be important in the non-sequential packs. 
One reason the non-continuum behavior described in Sec- 
tion |IV] is possible is that arbitrarily large forces are al- 
lowed; we can have a large positive force since there can 
be a negative force of a similar size to balance it. If 
no tensile contacts are allowed, then individual contact 
forces are limited in size by the magnitude of the applied 
external force. Unfortunately, because there is a unique 
solution to Equation [9] it is not possible to restrict the 
sign without making some fundamental changes to the 
packing geometry. Such changes did not affect the light 
cones in the sequential packs, but they may affect the 
results in the systems which did not average. Creating 
isostatic packs with tensile forces in a laboratory setting 
is also difficult. Tensile forces imply an attractive poten- 
tial, which will tend to cause particles to clump together 
and produce more contacts than would be present in an 
isostatic pack. An additional experimental complication 
is that the arbitrarily large contact forces which are al- 
lowed can drive a rearrangement of the system. Indeed, 
even if there are no tensile contacts, the system will still 
be quite delicate, and small applied forces will tend to 
lead to small rearrangements. As the applied forces get 
larger, this can qualitatively change the distribution of 
contact forces [50] . 

Despite the difficulties in observing these effects in a 
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laboratory setting, understanding the behavior of these 
systems is still important for the analysis of packs which 
are more readily observed. Even a system which is hy- 
perstatic can display isostatic behavior on small enough 
length scales. Furthermore, any transition from a dilute 
system to a densely packed one must pass through the 
isostatic state. 

Further work should be done to pinpoint what condi- 
tions lead to the existence of null stress behavior. Both 
an anisotropic contact angle distribution and isostatic- 
ity seem to be necessary but insufficient, as evidenced by 
the experiments in Section |VI| where the contact angles 
were changed while retaining the isostatic nature of the 
packs. We would like to improve our understanding of 
this issue by finding a better way to create non-lattice 
configurations with different contact angle distributions 
than the rather intrusive method described. It is possible 
that any sort of anisotropy in the pack construction will 
yield light cones. In this case, it would be interesting to 
to examine the connection between the light cones and 
network anisotropies to discover how the crossover from 
isotropic to anisotropic behavior occurs. 

A surprising result is that boundary modifications can 
evidently change the bulk force propagation from non- 
continuum to elastic, even at length scales smaller than 
the sample size, and even though only a vanishingly small 
fraction of the vibrational modes are affected. There 
must be a way to understand the contrasting force re- 
sponses in terms of these modes. We have begun prelim- 
inary investigations into the structure of the modes which 
are formed when an isostatic system is disconnected from 
its boundaries. These are a promising lead toward an 
explanation of both the presence of light cones and the 
dependence on the boundaries. 

Appendix A: Proof that V = MM T 

In Section |nj we claimed that the dynamical matrix T> 
and the contact matrix M are related via V = MM T , 
and offered a proof that depended on the invertibility of 
M. However, this relationship is true even if M is not 
invertible, as will be shown here. 

In the following, Greek indices will be used to denote 
contacts. Other indices will be written as a pair giving 
first the bead number and next the vector component. 
For example, for the displacement vector u, Ui X gives the 
x component of the displacement of the ith bead from 
its initial position. We also define cy to be 1 if beads i 
and j are in contact, and otherwise. c(i) is a list of the 
beads in contact with i; that is, all j such that c^- = 1. 

To begin, we will look at the structure of MM T , by 
computing (MM T \ x ,jx when i ^ j. This is just 

(MM T ) IXJX = ^M ix , a M jx>a 

From Scction|Hj we know that Mi X>a = cos 9ik if a is the 
contact between beads i and k, and otherwise. Thus 
Mix a Mj X a is only nonzero if i and j are in contact, 
and a is the contact between them. Each pair of beads 
can have at most one contact between them, so the sum 



collapses and we get 

(.)/.!/ 7 i, .......... = COS 9ij COS OjiCij = r<>> 2 (/,,<■, , 

If we were looking at the y components instead, we would 
get — sin 2 9ij Cij , and if we were looking at a mixture of 
x and y components, it would be — sin^- cos 0ycy. 

If i = j, then ^ Q M ixa M ixa will have a nonzero term 
for every contact that % has. We get 

feec(i) 

As before, if we were looking at y components we would 
have a sin 2 9^ term instead, and mixed components 
would give sin 9ik cos 9^. 

Now we want to compare this to the form of T>. From 
Equation 1 10| we can expand to sec 

C ~ y ^ \U'i X V>j X L'i X J X "I" HiylijyUyyjy 

id 

~\~ ^i X ^jy^i X: jy ~\~ ^iy^jx^iy.jx^ 

We split this into diagonal and off-diagonal parts, and 
modify the off-diagonal portion so that every pair is 
only written once: 

£ — ^ ^ (2Ui X Uj X T)i X j X -j- 'ZUiyTijyDiy jy (^^) 

i<j 

~t~ ^^ix^jy^ix.jy "T" '^V'iyV'j x 1~^iy,j x ^ 
k \ 2 2 

i 

Here we have used the fact that T> is symmetric, mean- 
ing *Dix,jy — ^jy^ix- 

Next we compute the energy between a pair of beads 
i and j when they are displaced some small amount. If 
they were not in contact before this displacement, then 
there is no interaction energy. If they were, then the 
energy is given by 

£ij = l(r ij -R ij ) 2 (A2) 

where Tlij is the vector pointing from bead i to j be- 
fore the displacements u 4 and Uj have been taken into 
account, and = + Uj — Ui is the same vector after 
the displacements have been made. The length of is 
given by 

rfj = (//,, (•(.>«,, • ,/,, a,, ■ ( //,, >\n(),, ■ a,,, //, v ) 

Taking the square root and assuming that the u.i are 
small compared to Rij , we see 

rij = Rij + cos 9ij {uj x -u ix ) + sii\9 VJ {u jv - u iy ) 

Putting this into Equation |A2[ we can sum over all 
pairs of beads to get the total energy, and expand it into 
the form 
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£ = - \2u lx u jx {- cos 2 OijCij) + 2u iv Ujy(- sin 2 %c^) (A3) 

i<j 

+ 2u ix u jy (- sin 6^ cos^cy) + 2u ly u jX {- sin% cosflyCy)) 
+ 2 22 + ^ COs2 + (*4 + u %) sin2 /y 'J r ',' 

+ 2(u tx u iy + u jx u jy ) s'mdij cos%Cjj J 



Now consider one term of the second sum, 



,, , u^)cos 2 c/,,( 



This is essentially a sum over contacts; for every 
pair in contact we get two terms with a factor of cos 2 % : 
one for bead i, and one for bead j. We can write this 
instead as a sum over beads, using the fact that cos 9ij = 
cos 2 9u- 



E 



The other two terms in the final sum of Equation [A3] can 
also be rewritten in this manner. Since this is true for all 
small vectors u, comparing the result with Equation |A1| 
gives us the terms of the dynamical matrix. We find 



V 
V 



"l-'-V = / y 



cos 2 9, 



sin 2 9 ik 



(A4) 



V — D 



sm f/ ifc cos « ifc 



feec(i) 



cos 2 OijCij 



^iytjy — sin 9ijCij 

These are exactly the same as the the elements of MM T 
that we found earlier. Thus even if M is not invertible, 
we still have V = MM T . 



Appendix B: Prevalence of boundaries which 
preserve isostaticity 

Here we give an argument to support the claim from 
the end of Section IIVI that even as the number of beads 
in a pack gets large, one can always expect to find a 
floor which does not create a contact matrix M with a 
null space. To do this, we estimate the probability that 
a given floor will produce no null vectors, and compare 
this to the total number of floors possible. 



Looking back at Figure [TTJ we can note that the 
nonzero elements of the null vectors are localized not only 
close to the floor, but also in different, possibly overlap- 
ping, horizontal areas. By deforming the floor in one 
place, we can eliminate the null vector located there and 
decrease the dimension of the null space without much 
affecting what happens along a different part of the floor. 
This indicates that it is some particular local configura- 
tions of floor beads which create the null vectors. We 
can model this by assuming that there is some probabil- 
ity per unit length a of a floor having a configuration 
which creates a null vector. 

Under this model, if many floors of a given length I 
are examined, the dimensions of their null spaces should 
follow a Poisson distribution with mean equal to a£. This 
is precisely what is shown in Figure [8] Using the fact 
that £ ~ vN, we find that a w 0.3 for all pack sizes. In 
units of bead-diameters, if the floor has length £, then 
the probability that there will be no null vectors is 



P(no null vectors) 







a°(l~a) e = (l-a) e (Bl) 



We next estimate the number of floors possible in a 
given pack. We limit ourselves to floors consisting of a 
chain of beads in contact with each other which wraps 
around the horizontal periodic boundary exactly once be- 
fore looping back on itself. Given some starting bead, we 
can approximate this as a one dimensional random walk 
where each step moves one bead to the right, and then 
randomly up or down. This is a crude approximation of 
the geometry in an isostatic pack, where each bead has 
on average four contacts; on either side, it will tend to 
have one contact higher in the pack, and one lower. This 
is exactly the situation for the diagonal isostatic lattice 



shown in Figure 19 1. 

After £ such steps, our path must return to the hori- 
zontal position of the starting bead, within one vertical 
step of the original height. In a one dimensional random 
walk, the number of paths of length £ which end at height 
k is 

N ^-{(£ + k)/2 
Using Stirling's formula, we find that the total number 
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of possible floors Np is 



i 

e/2 



(£-l)/2 



Combining this with Equation |B1[ we can compute the 
expected number of floors with no null space Eq to be 

E ~ (1 - aftr 1 '** 

Since a < 1/2, this grows with increasing path length. 
The path length i must grow at least as fast as the linear 
size of the system, so we can take I ~ y/~N. Thus the 
expected number of floors that do not create a null space 
in M grows with increasing system size. 



Appendix C: Stress fields in elastic media 

Here we state the results of Leonforte, et. al. [38] 
which we used in Section W\ to compare the stresses in 



isotropic packs with added hard ceilings with the stresses 
in an elastic solid. 

We will use an L x L box centered at the origin which 
is periodic in the x direction. The method is to solve 
the equations for linear elasticity assuming perfectly rigid 
hard ceilings and bottoms that are very rough, so that 
the displacements u(x,±L/2) there are zero. The ex- 
ternal force is modeled as a vertical external pressure 
p(x) in the middle of the pack. The medium is di- 
vided into two parts: part 1, for y > and part 2, for 
y < 0. The additional boundary conditions are then 

u«(x,0) = u( 2 )(x,0), o$(z,0) = o$-p(x), and 
Uxy{x, 0) = a£y(x,0). The finite boundaries quantize 
frequencies to multiples of Aq — (2%/L), and we can 
write the pressure as 



p(x) = E COs(q n x)s(q n )Aq 



(CI) 



71=0 



where q n = nAq. 

Instead of using a point force, we will choose p to be 
a Gaussian that is narrow enough that almost the entire 
area is contained in a length of approximately the same 
size as a bead diameter. Modest changes to this width 
do not much affect the results. 

Solving for the stress gives 



+ 4 V 2) = E C0S (<^) [« (1 ' 2) e 9 " a + b^e-M] 

oo 

^' 2) " ° { y V 2) = E cos(q n x)(q n y[a^e^ - b^e~^] + 2[c< 1 ' 2 V»* ~ d^e~^]) 

71=0 
OO 

°i» a) = E sin (<^) (^f t^ 1 ' 2 ^ 9 "" + 6 (1 ' 2) e" 9 " y ] + [cMe*** + d^e~^] ) 
'[(l + v)gL + 3- v] e- qL + {v - 3)e' 2<!L ] 
[(l + v)qL-3 + v]e- qL - (iz-3)^ 



a\q) 
a\q) 

b 2 (q) 
c\q) 

c 2 (q) 



n=0 
= s(q)Aq 
4S(q) 
= s(q)Aq 
4%) 

= w([- (1 + ^ L + 3 -^ 6 " L + ( ^ 3) ) 

= 8(l+lt%) ( [(1 + " )VL2/2 + (1 " V2)qL + (1 ~ V){3 " 
+ (1 - - 3)e~ 2qL ^j =d 2 (q) 

- sffT^S) (t (1 + " )2 ' 2i2/2 - (1 - - 2 >" L + (1 - "> (3 - 

+ (!-„)(„ -3)) =<!'(,)) 
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